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►^ . The Cherenkov radio pulse emitted by hadronic showers of energies in the EeV 

^S| I range in ice is calculated for the first time using full three dimensional simulations 

jn ■ of both shower development and the coherent radio pulse emitted as the excess 

charge develops in the shower. A Monte Carlo, ZHAire S, has been developed for 

this purpose combining the high energy hadronic interaction capabilities of AIRES, 

^ ■ and the dense media propagation capabilities of TIERRAS, with the precise low 

^D ■ energy tracking and specific algorithms developed to calculate the radio emission in 

ZHS. A thinning technique is implemented to allow the simulation of radio pulses 

induced by showers up to 10 EeV in ice. The code is validated comparing the results 

^ ■ for electromagnetic and hadronic showers to those obtained with GEANT4 and ZHS 

H . codes. The contribution to the pulse of other shower particles in addition to electrons 

and positrons, mainly protons, pions and muons, is found to be below 3 % for 10 

PeV and above proton induced showers. The characteristics of hadronic showers 

and the corresponding Cherenkov frequency spectra are compared with those from 

purely electromagnetic showers. The dependence of the spectra on shower energy 

and high-energy hadronic model is addressed and parameterizations for the radio 

emission in hadronic showers in ice are given for practical applications. 
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1 Introduction 



The observation of neutrinos of EeV-scale energies is one of the main priorities 
in Astroparticle Physics. The detection of neutrinos will open a new window to 
observe parts of the Universe shielded by large depths of matter, unaccessible 
to conventional astronomy. The measurement of diffuse neutrino fluxes will 
provide further clues to the identification of the sources of ultra-high energy 
cosmic ray production, their composition and their production mechanisms pLj. 
In addition such detections could have important implications for fundamental 
particle physics. 

A very promising and cost-effective method to detect high-energy neutrino 
interactions was first proposed by G. A. Askar'yan in the 1960's p]. The idea 
is to detect the Cherenkov radiation at radio wavelengths generated by the 
excess number of electrons in the cascade of particles resulting from a high- 
energy particle interaction in a dense medium transparent to radio waves. 
The development of this excess charge, due to the interactions with matter 
electrons, is often referred to as the Askar'yan effect. The effect has been 
experimentally confirmed in accelerator experiments at SLAG in media such 
as sand [3'3], rock salt [5] and ice [6], with results in good agreement with 
theoretical calculations [7]. At radio frequencies (MHz-GHz), the emission is 
coherent and the radiated power scales with the square of the primary particle 
energy. This makes this method very promising for the detection of neutrinos 
and cosmic rays of the highest energies (EeV) PI^PH] . 

Several experiments have already exploited this technique searching for ultra- 
high energy neutrinos, but no positive detection has been reported so far. 
These include experiments using the ice cap at the South Pole, such as the 
ANITA balloon experiment [TT|12fl3] and the RIGE array of antennas buried 
under the Antarctic ice [T^fT^ . Other arrays are beginning to be constructed, 
such as the Askar'yan Radio Array ARA [T6] and ARIANNA [17]. There are 
also a number of projects using the Moon as target for neutrino interactions, 
along with radio telescopes on Earth as radiation detectors, such as the pio- 
neering Parkes [TM^ . GLUE [20], Kalyazin [51], NuMoon [22], LUNASKA 
and RESUN 



The interpretation of data from these experiments requires a detailed knowl- 
edge of the magnitude, angular distribution and frequency-dependence of the 
emitted Gherenkov radiation, which can then be related back to the proper- 
ties of the induced cascade. This calls for accurate simulations of the Fourier- 
spectrum of coherent Gherenkov radiation from EeV showers in different dense 
media. In the past, full simulations of electromagnetic (EM) showers in ice 
were done up to PeV energies, using the well-known and well-tested ZHS code 
[5)125] and different versions of the GEANT code [25P7] . Also, full simulations 



of hadronic and neutrino-induced showers in ice were carried out up to 100 
TeV with GEANT [28]. Hybrid Monte Carlo [29] and thinning techniques 
were also developed to simulate electromagnetic [ST|32f34] . hadronic 



and neutrino-induced showers [M|55] above PeV energies (mainly in ice and 
water). Above these energies the LPM effect [36)137] starts to be important in 
the media under experimental consideration [37|38|31] . Semi-analytical calcu- 
lations of the radio-emission have also been performed [39|^ . Very recently, 
simulations of not only the Fourier-components of the spectrum but also of 
the radio-pulse in the time-domain emitted in electromagnetic showers, have 
been performed with the ZHS code[42j. 

Modeling hadronic showers up to EeV energies is of utmost importance for 
neutrino detection, since they are induced by all neutrino flavors in neutral 
current (NC) interactions, as well as in the hadronic vertex of charged current 
(CC) interactions of muon and tau neutrinos. Although at EeV energies only 
~ 20% of the neutrino energy is on average carried by the hadronic shower, 
these showers are known to be less affected by the LPM effect [33], and their 
Cherenkov emission does not suffer from the shrinking of the Cherenkov cone, 
which would otherwise reduce the solid angle for observation. On the other 
hand, mixed showers induced in CC electron neutrino interactions are com- 
posed of a purely electromagnetic shower, produced by an energetic electron 
carrying on average ~ 80% of the energy of the neutrino, and a hadronic 
shower, initiated by the debris of the interacting nucleon. The advantage for 
neutrino detection of these types of interaction is that all the neutrino en- 
ergy is channeled into the resulting shower. However, the observation of the 
electromagnetic shower is expected to be very difficult, since the LPM effect 
shrinks the angular distribution of the electric field, reducing dramatically the 
available solid angle for detection. For these reasons, experiments aiming at 
detecting neutrinos using the radio technique gain most of their acceptance 
from neutrino-induced hadronic showers 



Clearly, accurate simulations of hadronic showers with energies above which 
the LPM effect becomes effective in different dense media are needed. In this 
paper we present the first steps in that direction. We have used the well-known 
AIRES code [H] in combination with the TIERRAS package [12] to simulate 
proton-induced showers in dense media, such as ice. We have implemented the 
algorithms to calculate the Fourier components of the electric field produced 
by charged particle tracks in the shower, developed by Zas-Halzen-Stanev in 
the well-known and well-tested ZHS code [9||25] . The result is a flexible and 
powerful code named ZHAire S, which allows the simulation of electromag- 
netic, hadronic and neutrino-induced showers in a variety of media, along with 
their associated coherent radio emission due to the Askar'yan effect. 

It is important to stress that in [3T|33] the calculations of the emission prop- 
erties of showers in ice up to EeV energies were only done in an approximate 



way. They were based on the Fourier-transform of simulated longitudinal and 
lateral shower profiles. With ZHAire S we can obtain the radio-pulse features 
in a consistent manner within a well-tested simulation, computing the electric 
field emitted by charged particles in the shower. In Ref . [28j . GEANT simu- 
lations were used to calculate the radio-pulse properties in the same way as 
the ZHAire S code. However these calculations were limited to energies below 
100 TeV due to GEANT limitations, while with ZHAire S we can simulate 
showers at EeV energies and above. Also, and for the first time, the contribu- 
tion to the radio pulse of charged pions, muons and protons is accounted for 
and quantified, and parameterizations of the frequency spectrum of the radio 
pulse due to Askar'yan effect in hadronic showers are presented. Previous pa- 
rameterizations of hadronic showers exist in the literature [33], but these are 
based on less accurate simulations than those presented in this work. 

This paper is structured as follows. In Section [2] we describe the new ZHAire S 
code and compare its performance and results for electromagnetic showers with 
those of the well-tested codes GEANT 4 p] and ZHS [251132] . We also explore 
thinning techniques in the ZHAire S code, essential for the simulation of EeV 
showers. Section [3] is devoted to the simulation of hadronic showers up to EeV 
energies, emphasizing their differences with respect to purely electromagnetic 
showers. Section H] summarizes and concludes the paper. In Appendix A, we 
also give parameterizations of the intensity of the Fourier-spectrum of the 
Cherenkov electric field emitted in hadronic showers, as a function of shower 
energy and observation angle. 

In this paper we concentrate on hadronic showers in ice. We defer to future 
papers the study of coherent radio emission from EM and neutrino-induced 
showers. However we show that our results for proton showers can be applied 
to approximately model the Cherenkov emission from hadronic showers in- 
duced in high-energy neutrino interactions. The code can also be applied to 
the simulation of radio emission in extensive air showers [50|5T] . despite the 
different emission mechanism 



Radio emission in electromagnetic showers using the ZHAire S 
code 



Coherent Cherenkov emission from high-energy showers in a dense medium 
such as ice is due to the excess of electrons over positrons in the shower. This 
excess is mainly caused by the low energy scattering processes, such as Comp- 
ton, M0ller and Bhabha scattering, that incorporate electrons of the medium 
into the shower [25], as well as to positron annihilation. The characteristics 
of this emission are dependent on the lateral and longitudinal profiles of the 
shower[7]. The longitudinal profile affects mainly the emission at angles far 



from the Cherenkov angle. As the shape of this profile is governed by high en- 
ergy interactions, radio emission is strongly affected by the LPM effect, which 
reduces the probability of high energy {E > 1 PeV in ice) EM interactions. 
The lateral profile affects the features of the spectrum close to the Cherenkov 
angle at frequencies typically above ~ 1 GHz, and is mainly determined by 
low energy particle scattering. 

Next we describe the key elements of the ZHAire S code, a flexible and pow- 
erful Monte Carlo that allows the simulation of electromagnetic, hadronic and 
neutrino-induced showers and their associated coherent Cherenkov radio emis- 
sion due to the Askar'yan effect in a variety of media. 



2.1 The ZHAire S code 

ZHAire S is based on the well-known AIREq^ code [Sj which we have used 
in combination with the TIERRAS package [15] to simulate showers in dense 
media, such as ice. We have implemented in AIRES algorithms to calculate the 
Fourier components of the electric field produced by charged particle tracks 
in the shower. These algorithms are the same used in the well-known and 
well-tested ZHS code, developed by Zas, Halzen and Stanev pi25] . 

The standard procedure to simulate Cherenkov coherent radio emission due 
to the Askar'yan effect in an energetic shower is to divide all charged particle 
tracks in the shower, down to a relatively low threshold energy, in small steps 
which are approximated by straight lines in which particles travel at constant 
speed. A space-time position is associated with the end points of each step, 
along with the speed and the charge of the particle. This is all that is needed 
to calculate the electric field, and to properly take into account interference 
effects between all steps (for a detailed discussion and the relevant equations 
see P?P7] ). The AIRES package [H] is a well established Monte Carlo ex- 
tensive air shower simulation code, from which detailed information on the 
space-time positions of the tracked particles can be extracted. This allows the 
use of AIRES (v2.8.4a) to simulate the coherent radio emission due to any type 
of particle interaction. However, AIRES by itself can only simulate showers in 
air, where the contribution of the Askar'yan mechanism to radio emission is 
thought to be small, and other processes, such as geo-synchrotron radiation, 
are expected to be more important [H]. In order to simulate showers in dense, 
dielectric media, we have used the TIERRAS package |l5j . Originally designed 
to continue the simulation of an extensive air shower underground, TIERRAS 
also enables the simulation of particle cascades in other dense media, such as 
ice, sea water, rock and soil. TIERRAS accounts for the atomic and mass num- 
ber dependence of the interaction cross-sections and energy-losses of baryons. 



AIRshower Extended Simulations 



mesons and leptons, and includes a detailed treatment of the LPM effect of 
special relevance for radio emission [38|31f33] . 

In this work we concentrate on showers developing in ice, with refraction index 
n = 1.78 and density p = 0.924 g cm'^. Using the ZHAire S + TIERRAS 
framework, we track e^ with energies above 80 keV, and tt"^, fi"^, p andp above 
500 keV, neglecting radio emission for energies below these thresholds [fj. These 
energy cuts are the lowest available in AIRES. We also do not account for the 
emission of other charged particles such as kaons, which as will be shown in 
Section [3] constitute a negligible contribution. 

We used the minimum propagation step available in AIRES, which is typically 
much smaller than a radiation length. Directly from the AIRES routines we 
get the start and end points for each step, along with the charge of the particle 
and propagation time. From these tracks the radio emission is then calculated 
making extensive use of the formula for the frequency (u) spectrum of the 
electric field E, at a position x in the Fraunhofer limit of large observation 
distance R, derived in PE5] and reproduced here for self-consistency: 



E(u;,x) = ^^^^e'^'^-'^-^^'^v^ 



27reoc2 R 



i{uj — k ■ v) 



where /Xr is the relative magnetic permeability of the medium, eo the free space 
permittivity, c the speed of light in vacuum, k = un/c is the wave vector in 
a medium with refractive index n, v is the particle velocity, ti the start time 
and St is the time interval between the end points of the track. 

We recall that this equation has been obtained with the following convention 
for the Fourier transform of the electric field: 



Ru^) 



f{t) e"^'dt 



(2) 



where the factor 2 corresponds to an unusual convention (this factor is usually 
either 1 or (27r)~2). This does not pose a problem as long as the inverse Fourier 
transform is done consistently (see [12] for more details). 

The necessary division of the particle paths into piecewise linear tracks for 
radio pulse calculations can be performed in various ways, with shorter divi- 
sions allowing increasingly accurate parameterizations to higher frequencies at 
the expense of computer time |27]. In ZHAire S the splitting is finer than in 



^ Note that due to the isotropically distributed directions of particles below these 
thresholds, their contribution to the electric field is expected to be very small [52] . 



the default ZHS [39], roughly 3-4 times more steps, which is accurate enough 
for our purposes up to ~10 GHz and even above, i.e. above the frequency at 
which the coherent behavior of the emitted signal is expected to break down 
in ice [7]. Shorter tracks will lead to a significant increase in computing time 
but at no improvement in accuracy. 

Coherent Cherenkov radio emission in showers in dense media is mainly due 
to electrons and positrons with ~MeV kinetic energies. For instance, 50% of 
the projected track length due to the excess of negative charge in the shower, 
a quantity known to be proportional to the coherent electric field, is produced 
by electrons with kinetic energy between 100 keV and ~ 6 MeV, with little 
contribution expected below 100 keV [23]. We have exphcitely checked this 
in our simulations, and the results obtained agree with [25]. A very accurate 
treatment of this keV to MeV kinetic energy range is thus needed to obtain a 
precise determination of the radio emission intensity. However, in AIRES, elec- 
trons from knock-on (KO) electron and positron interactions, such as Bhabha 
and M0ller scattering, are not tracked explicitly if the kinetic energy of the 
incident e^ is below Kiow = 1 — 2 MeV, and a continuous energy loss dE/dX 
is used instead. This posed a difficulty, since the approximation was not ac- 
curate enough to describe the sub-MeV electrons, and also because in order 
to fully describe the low energy part of the shower electronic component, 
one should lower -ft'iow to explicitly simulate knock-on interactions for elec- 
trons down to -ft^iow = 106 keV. These two issues made the parameterization 
of dE/dX used in AIRES unsuitable for radio applications. To solve these 
problems, we calculate the continuous ionization loss explicitly for low energy 
electrons and positrons integrating the energy loss for secondaries of energy 
below Kiow = 106 keV, as done in the ZHS code, and replace the dE/dX 
parametrization implemented in AIRES with the actual Bethe-Heitler energy 
loss formula at low energies. We have checked that the dE/dX obtained in 
this way is in very good agreement with that implemented in GEANT 4 
with maximum relative differences ~ 1% when the same -ft'iow is used. 



2.2 Comparison to GEANT4 and ZHS 



In order to validate the ZHAire S code we have performed several compar- 
isons of its output in ice with that produced by the GEANT 4 [3^ and ZHS 
[25] codes. We have concentrated our efforts in comparing the lateral (per- 
pendicular to shower axis) development of the excess negative charge in the 
shower, as well as the track-length due to this excess negative charge. For 
these comparisons, we raised the e^ kinectic energy cuts to 106 keV, to match 
the thresholds used in GEANT4 and ZHS. 

Two relevant quantities affecting significantly the longitudinal development of 
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Fig. 1. Top panel: Mean free path (A) of photons at different energies in ZHS (long 
dashes) and ZHAire S with (short dashes) and without (solid line) photoproduc- 
tion interactions enabled. Also shown is the interaction length in CORSIKA [3l] 
with both pair production and photoproduction enabled (triangles). Note that pho- 
toproduction is not accounted for in ZHS. Bottom panel: Mean free path of electrons 
of different energies in ZHS (blue solid line) and ZHAire S (red dashed line over- 
lapping the blue solid line). 



an electromagnetic shower are the pair production (PP) and bremsstrahlung 
(BR) cross sections. In order to compare the PP cross-sections in ZHS and 
ZHAiRE j'^l. we obtained the photon mean free path (MFP) at several ener- 
gies. The result is shown in the top panel of Fig. [TJ Since photoproduction is 
not taken into account in ZHS, we disabled it temporarily in some ZHAire S 
simulations in order to allow a direct comparison. Up to 1 PeV, the MFP for 
PP interactions obtained in ZHAire S is very similar to that in ZHS, with dif- 
ferences ~ 3%, and is equal to (9/7)Xo ~ 46 g cm~^, where Xq = 36.08 g cm~^ 
is the radiation length in ice. Above ~ 1 PeV, the energy at which the LPM 
effect starts to produce a decrease in the PP cross section in ice, the aver- 
age interaction depths begin to deviate from each other. For photon energies 
above 100 PeV, the ZHAire S PP cross-section is around 30% lower than the 
ZHS one, and the LPM effect turns on at a smaller energy in the ZHAire S 
(TIERRAS) code than in ZHS. We attribute this difference to details in the 
implementation of the LPM effect in both codes [551l25f37] . This will not be 
investigated further in this paper, since these differences are mostly irrelevant 
for the hadronic showers and radio pulses discussed below |33]. Above ener- 
gies of order 1 EeV, the photoproduction cross-section becomes larger than 
the PP cross-section and the interaction length decreases with respect to the 
case in which only the LPM PP cross-section was accounted for. The energy 
at which this happens is however model-dependent [551156] . The average in- 
teraction length when photoproduction interactions are enabled in ZHAire S, 
and in CORSIKA are also shown for comparison. 

The difference in the PP cross section is partly responsible for differences in 
some macroscopic shower observables, such as the position X^ax at which the 
number of particles in the shower is maximum. In Fig. [2] we show the average 
Xmax in electron-initiated showers simulated with ZHS and ZHAire S (with 
and without photoproduction interactions enabled). The PP cross-section at 
high energies is lower in ZHAire S due to the difference in the strength of 
the LPM effect. This translates into an up to ~ 20% larger average X^ax 
in ZHAire S. If we disable the LPM effect and photoreactions in both ZHS 
and ZHAire S, the elongation rate SXmax/f^ logio -E becomes a constant, as 
expected, and we obtained a value of 81 g cm~^ in both cases. 

In the bottom panel of Fig. [T] we also compare the mean free path for bremsstrahlung 
of electrons in ZHAire S and ZHS simulations which agree to better than 

-5%. 

In Fig. [3] we show the lateral distribution of the excess negative charge, at a 
depth around shower maximum, obtained from simulations of 1 PeV electron- 
induced showers in ice, performed using ZHAire S and ZHS. The lateral dis- 



^ Note that ZHAiRE S uses the same high energy cross-sections as the 
AIRES +TIERRAS framework, including the treatment of the LPM effect. 
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Fig. 2. Average depth of shower maximum Xmax in electron-initiated showers simu- 
lated with ZHS (dark red dots) and ZHAire S with (light green dots) and without 
(asterisks) photoproduction enabled. The dashed lines, shown just as a reference, 
are fits to ZHS data: a power law at high energy and a logarithmic lengthening 
(elongation rate) at lower energies. 

tribution is a relevant quantity that mainly affects the frequency spectrum of 
the field at high frequencies, typically above ~ 1 GHz. The agreement between 
both codes is rather gooco 

In Table [1] we show the average track lengths obtained from 100 simulations 
of 100 GeV electron-induced showers with ZHAire S and ZHS, and compare 
them with GEANT results taken from [27J. The total track length is defined as 
the sum of all electron and positron tracks in the shower. The total projected 
track length is the sum of all the projections of these tracks onto the shower 
axis. The excess projected track length is the difference between the sum of 
all electron track projections minus the sum of all positron track projections. 
This is the most relevant quantity that determines the normalization of the 
electric field spectrum, since electrons and positrons contribute to the field 
with opposite signs. The relative differences between the total track length 
and excess projected track length as obtained in ZHAire S and ZHS are 4.1% 
and 0.9%, respectively. The correspondin g d ifferences between ZHAire S and 
GEANT4 are 3.6% and 1.0%, respectiveljE]. 



^ The differences near shower axis are mostly due to the fact that ZHS by default 

includes all particles with r < 0.1 g cm~^ in the first bin. 

^ Relative differences are calculated using [ZHAire S - ZHS] /ZHS and [ZHAire S 

- GEANT4]/GEANT4. 
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Fig. 3. Distribution of the excess negative charge as a function of the distance to 
the shower axis, obtained from simulations of 1 PeV electron-induced showers in ice 
using ZHS (red dashed line) and ZHAire S (blue solid line) . The lateral distribution 
refers to a depth around Xmax ~ 600 g cm~^. See text for discussion. 
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MC Code 



ZHAire S GEANT 3.21 GEANT 4 ZHS 



Total Track [m] 566.8 

Total Projected [m] 459.3 

Excess Projected [m] 123.9 

Excess Projected/Total 0.219 



577.9 



450.0 
123.5 

0.214 



587.9 



453.2 
122.7 
0.209 



591.1 



480.9 
125.1 
0.212 



Table 1 

Comparison between average track lengths obtained from simulations of electron- 
induced showers of ^o =100 GeV in ice, using ZHAire S, GEANT 3.21, GEANT 4 
and the latest version of ZHS [32] . The results of ZHAire S and ZHS are the average 
of 100 simulated showers each, while the results of the GEANT 3.21 and GEANT 
4 codes are taken from [2^. The total excess track length is 193.4 m in ZHAire S 
and 200.3 m in ZHS. Corresponding total excess track lengths for GEANT are not 
available. Note that the results shown here for the ZHS code give better agree- 
ment than reported in [27] . due to a minor correction in the bremsstrahlung terms 
proportional to the density of matter electrons. 

In Fig. m we compare the average frequency spectra of the Cherenkov radiation 
obtained from 20 ZHS simulations of electron-induced showers of primary 
energy Eq = \ PeV, with the corresponding average spectrum obtained with 
the ZHAire S code. One can see a very good agreement below the cut-off 
frequency in the spectrum (the frequency at which the electric field peaks 
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Fig. 4. Average frequency spectrum of the Cherenkov radiation obtained in simu- 
lations of 20 electron-induced showers with primary energy Eq = 1 PeV, using ZHS 
(blue solid line) and ZHAire S (red dashed line). The spectrum is shown at three 
observation angles with respect to the shower axis. The RMS of the 20 simulated 
showers is also shown. 



before decreasing). Beyond the cut-off frequency the agreement worsens. At 
the Cherenkov angle and frequencies above about 3 GHz, the electric field is 
sensitive to the fine details of the shower, in particular to the different splitting 
of the charged particle tracks in ZHS and ZHAire S codes. The finer splitting 
in the ZHAire S code (3-4 times more steps) compared to ZHS [32] in which 
the electric field is calculated within the approximation "b" [39], leads to a 
higher spectrum in the GHz frequency range, when compared to ZHS. The 
effect of the splitting on the high-frequency Fourier components of the electric 
field was demonstrated numerically in ^39j. At higher energies, typically above 
1 EeV or so photoproduction interactions start to play an important role in 
ice (see Fig. [2] and Ref. [lO])- While in the ZHS code only EM processes are 
accounted for we expect the results of the ZHAireS and ZHS code for radio 
emission in EM showers to differ. This topic will be studied in a forthcoming 
paper [52] . 
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2.3 Thinning in ZHAire S 



In order to make possible the simulation of showers above PeV energies, AIRES 
uses thinning algorithms [30lH^ to reduce computing time significantly. The 
basic idea of thinning is to follow only a small, representative fraction of the 
particles in a shower, and assign to each tracked particle a corresponding 
weight to compensate for the rejected ones |30j . 



Based on the properties of radio emission, a study was made in [32| of a 
thinning technique for calculations of radio emission in EM showers. The con- 
clusion was that applying a thinning algorithm similar to the original Hillas 
algorithm [30], but with two thinning thresholds instead of one, i.e., thinning 
only particles with energies between -Emin and E^nx, produces an accurate 
representation of the coherent Cherenkov radio-emission in electromagnetic 
showers, while at the same time keeps computing time to an acceptable value. 
The thinning parameters obtained in |32] were E'min ~ 100 MeV — 1 GeV 
and -Emax ~ (10~^ — 10~^) Eq, where Eq is shower energy. These parameters 
were obtained as a compromise between accuracy (smaller than ~ 10%) and 
CPU time. This algorithm is now implemented in the ZHS code. 

The ZHAire S code uses the same thinning algorithm as AIRES, with an en- 
ergy threshold E'th, equivalent to i?max in ZHS, and a so-called weight factor 
Wf [H], which sets a maximum weight allowed in the simulation. It is im- 
portant to remark that the weight factor is not the maximum weight of the 
particles in the simulation, and that the maximum weight increases linearly 
with shower energy [H]. By fine-tuning Eth and Wj, one can make Wf work 
in a similar fashion as -Emin, thus making the ZHAire S thinning algorithm 
work in a similar way as the one in ZHS. 

The key point is to select thinning parameters {Ef^ and Wf, in the case of the 
ZHAire S code), so that the computing time is minimized, while maximiz- 
ing the agreement between full simulations, i.e., those performed following all 
particles, and the thinned ones. For that purpose we have compared full simu- 
lations of 20 showers at 1 PeV with thinned simulations using different choices 
of E'th and Wf, in a similar fashion to what was done in [32] • A good compro- 
mise between accuracy and computing time was obtained for the parameters 
Eth = 10""^ Eq (as in the case of ZHS [32j) and Wf = 0.06. An example 
is shown in Fig. [5l where we plot a comparison between 20 thinned simula- 
tions of electron-induced showers of energy 1 PeV performed with the thinning 
parameters above, and 20 full simulations of the same energy, both using the 
ZHAire S code. One can see that the agreement is very good up to frequencies 
beyond the cutoff frequency in the spectra, which depends on the observation 
angle. In this region, the difference between the average electric field in thinned 
and full simulations is always below 4%, and the difference in the RMS of the 
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Fig. 5. Average frequency spectrum of the Cherenkov radiation generated by elec- 
tron-induced showers with primary energy Eq = 1 PeV, obtained from 20 un-thinned 
simulations (red dashed line) and 20 thinned simulations (blue solid line) , using the 
ZHAiRE S code. The spectrum is shown at three observation angles with respect to 
shower axis. The RMS of the 20 simulated showers is also shown. The parameters 
used in the thinned simulation were £'th = 10~^ Eq and Wf = 0.06. 

simulated showers normalized to the average value, i.e., (cr^™ — a^^^^^) / {\E\) 
with {\E\) the electric field, is always smaller than ~ 5%. However, for obser- 
vation angles away from the Cherenkov angle and frequencies well beyond the 
cutoff in the spectrum, the two calculations deviate significantly. This is not 
surprising since the Fourier components are sensitive to the fine structure of 
the shower (even at the individual particle level) in this angular and frequency 
region, and clearly the thinning algorithm in the ZHAire S code is unable to 
account for it. This is not a problem because in that region the electric field 
strength is typically a factor 10 to 100 smaller than in the coherent region, 
and it is not expected to contribute significantly to the total emitted power. 

At energies above ~ 1 — 10 PeV, simulating showers without thinning becomes 
impractical, since computing time scales linearly with shower energy. This 
makes a direct comparison between thinned and un-thinned simulations above 
1-10 PeV almost impossible. Instead, we have compared ZHAire S thinned 
simulations, using the same parameters described before, with thinned sim- 
ulations performed with ZHS, using its own thinning parameters {E^i^ and 
-E'max); for which a comprehensive study of their validity was made [32]. The 
agreement between the two simulations is excellent at all energies, and we 
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Fig. 6. Same as Fig. [5] for proton-induced showers. 

adopt the thinning parameters (-Eth = 10~^ Eq, Wf = 0.06) in our subsequent 
simulations. 



3 Radio emission in hadronic showers 



In this section we study the coherent Cherenkov radio emission of hadronic 
showers in ice. Using the ZHAire S code, we simulated large samples of proton- 
induced showers up to an energy oi Eq =10 EeV. 

A first important check to perform is whether the values for the thinning pa- 
rameters ii^th and Wf, obtained using EM showers in section [2^3} are also valid 
to accurately calculate radio pulses from hadronic showers. For this purpose, 
in Fig. |6] we compare 20 full simulations of proton-induced showers with 20 
showers simulated using the same thinning parameters obtained for electro- 
magnetic showers. The agreement between both simulations is equivalently 
good except again for the frequency region above the region of full coherence. 
The CPU time required remains within reasonable values, reinforcing the con- 
clusion that the same thinning algorithms and parameters can be used safely 
and provide a good description of both, electromagnetic and hadronic showers. 
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3.1 Comparison to GEANT4 simulations 



In this section we compare our results to simulations of proton-induced show- 
ers using GEANT4, as reported in [28]. In Fig. [TJ we compare the projected 
excess track length obtained in simulations of proton-induced showers using 
GEANT4 and ZHAire S. For clarity, the projected excess track length has 
been divided by shower energy Eq. For the purposes of the discussion be- 
low, we also show in Fig. [7] the excess projected track length in the case of 
electron- induced showers, obtained in ZHAire S. 

GEANT4 simulations of proton showers have only been performed up to 90 
TeV, while we performed ZHAire S simulations up to 10 EeV, an energy 
5 orders of magnitude higher. One can see that there is a good agreement 
between the results of both simulations. The small differences between the 
track lengths from both simulations, up to a maximum of ~ 7%, are well within 
1 RMS of the ZHAire S results, which were obtained from 20 simulations for 
each primary energy. In the case of the GEANT4 simulations reported in |i28j, 
a RMS < 10% is quoted, except for 90 TeV, for which only a single shower 
was simulated and thus the RMS could not be obtained. 

A feature that can be seen in Fig. [7] is that the track length due to the excess 
of electrons in hadronic showers approaches that in electron-induced showers 
at high energies. This is due to the fact that the fraction of primary energy 
going into the electromagnetic component of the shower increases with shower 
energy. Assuming energy equipartition between pions, in each interaction a 
fraction of 1/3 of the incident pion energy would go into the electromagnetic 
component. Because the medium is dense, the produced charged pions are 
expected to interact before decaying and hence produce more vr^s that con- 
tribute further to the EM content of the shower. As shower energy increases 
more particle generations and interactions occur, and hence more energy is 
transferred to the EM component. Hadronic showers at the highest energies 
will have a markedly electromagnetic character because of the high energy in- 
volved and the medium density which prevents decays of most of the charged 
pions. 

To further check this interpretation we have explicitly obtained in ZHAire S 
simulations the energy going into the EM component of the shower. The result 
of this study shows that the increase with shower energy of the EM energy 
fraction is responsible for the increase with energy of the excess electron track 
in proton-induced showers shown in Fig. [71 Between 1 TeV and 90 TeV, the 
ratio between the GEANT4 e^ track length calculation, as shown in figure 
[TJ and the EM energy fraction obtained from ZHAire S simulations varies 
only slightly, between 11.6 and 11.8 m/TeV, showing that the projected e^ 
excess track length is proportional to the EM energy fraction of the shower. 
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Fig. 7. Projected excess electron track length as a function of energy, obtained in 
GEANT4 [2H] (solid line) and ZHAire S (dot-dashed Hne) proton-induced shower 
simulations. For clarity, the track length has been divided by shower energy. Also 
shown are the track length for electron-induced showers, simulated using ZHAire S. 
The GEANT4 results have an uncertainty (not plotted) of ~ 10% [28j. 

We have also checked that the total track length behaves similarly, i.e. it is 
proportional to the EM energy fraction, confirming that the behavior seen in 
figure [7] is not an effect of projecting the track along the shower axis. The EM 
energy fractions obtained here are also in agreement with those obtained in 
[33], where a simplified simulation of shower development was performed. 

Interestingly, in the case of electron-induced showers, one can see in Fig. [7] a 
small deviation from the linear dependence of the track length with energy. 
We attribute this behavior to photonuclear interactions becoming relevant at 
high energy and contributing to the hadronic component (i.e. reducing the 
EM energy fraction), even in electron- induced showers. This is confirmed by 
our simulations, as can be seen by the values of the EM energy fraction in 
electron induced showers, shown in table [2j 



3.2 Radio emission 



The dominant radiation mechanism in hadronic showers is the same as in 
electron-induced ones, namely, emission of coherent Cherenkov radiation from 
the excess of electrons over positrons, i.e., the Askar'yan effect. However, there 
are differences between hadronic and electromagnetic showers in the normal- 
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ization of the magnitude of the electric field, as well as in some important 
features of the frequency spectrum, which stem mainly from the different spa- 
tial distribution and from the fraction of EM energy in both types of showers 



FigJT] displays the projected track length excess divided by primary shower 
energy as a function of shower energy for both electromagnetic and hadronic 
showers. This ratio is almost constant for electromagnetic showers while for 
hadronic showers it increases with energy because the EM energy ratio of 
hadronic showers also increases with energy, as explained before, approach- 
ing the value for electromagnetic showers but remaining always smaller. This 
reflects itself in the normalization of the electric field which is smaller in 
hadronic showers. This can be appreciated in the top panel of Fig El where 
the Cherenkov spectrum of proton and electromagnetic 1 PeV showers are 
compared. At this energy, the ratio between the electric fields at 1 MHz and 
at the Cherenkov angle of the proton and electron induced showers is 0.864, 
whilst the ratio between the EM energy fraction of these showers is very sim- 
ilar, 0.902, i.e. {\E\p) / {\E\e) ~ ^m/™^*^°"/EM/™'=*^°". It can be also noted 
that as the shower energy increases (from top to bottom panels), the electric 
field produced by proton-induced showers approaches that of electromagnetic 
showers at small frequencies, where the coherence is full (10 - 100 MHz de- 
pending on the observation angle), as expected from the results of the excess 
projected track length displayed in Fig [71 

A relevant feature of the frequency spectrum is the cut-off frequency at the 
Cherenkov angle, which is known to be inversely proportional to the lateral 
dimensions of the shower |7]. Hadronic showers typically spread over a larger 
lateral distance than electromagnetic ones [28], and as a consequence, the cut- 
off frequency in the hadronic case is expected to be smaller. The Cherenkov 
spectrum of electron-induced showers peaks at ~ 3 GHz for observation at 
the Cherenkov angle, while the peak shifts to lower frequencies, around 2.5 
GHz, in proton-induced showers. The effect is however difficult to appreciate 
in Fig m 

Another important effect is the shift in the cut-off frequency at angles away 
from the Cherenkov angle, which is known to be determined by the longitudi- 
nal development of the shower [7f25] . Electron- induced showers develop more 
slowly than hadronic ones and hence penetrate more in the medium. As a 
result, for a fixed energy, hadronic showers are on average shorter in the longi- 
tudinal direction (along the shower axis) and we expect the cut-off frequency 
of their spectra, away from the Cherenkov angle, to be typically larger than in 
electron-induced showers. This effect is however not visible at energies below 
the LPM scale (~ PeV) (top panel of Fig. |8]), but it is apparent at energies 
above which the LPM effect becomes important for shower development. 



Hadronic showers are known to be less affected by the LPM effect than electro- 
magnetic ones [33], the reason is that the production of high energy photons 
through neutral pion decays is suppressed above PeV energies, since these 
energetic pions are more likely to interact than to decay. As a result, while 
electromagnetic showers are dramatically stretched in the longitudinal direc- 
tion, hadronic showers grow much more slowly with energy. The immediate 
consequence is that, for a fixed observation angle away from the Cherenkov 
cone, the cut-off frequency is significantly smaller in electromagnetic than in 
hadronic showers. This is apparent in Fig. [SI At 100 PeV and 10° inside the 
Cherenkov cone, the cut-off frequency is ~ 50 MHz in the electron-induced 
showers, while it is ~ 80 MHz in the proton-induced ones. When EeV energies 
are reached, this effect becomes very strong. Moreover, shower-to- shower fluc- 
tuations of the longitudinal spread become very large due to the LPM effect 
in electron-induced showers |2Z]. As a result, some individual electromagnetic 
showers have a cut-off frequency even 10 times smaller than those in hadronic 
showers. In Fig. [SI the vertical bars represent the RMS of the average of 20 
showers, making apparent that the RMS is much larger in electromagnetic 
than in hadronic showers for angles away from the Cherenkov one, for essen- 
tially all frequencies. It can also be appreciated that fluctuations significantly 
increase as shower energy increases in the case of electron-induced showers, 
while in hadronic showers the dependence on energy is weak. 

We compare hadronic showers among themselves in Fig. [9], where we have 
divided the electric field by shower energy to make the differences between 
the features of the spectrum more apparent at different energies. The cut- 
off frequency at the Cherenkov angle increases slowly as the shower energy 
rises, since it is determined by the overall lateral spread of the shower, which 
becomes slightly narrower as the energy increases. Away from the Cherenkov 
angle, one can clearly see that the cut-off frequency decreases with energy as 
expected, but not so rapidly as in electromagnetic showers, because in hadronic 
showers the longitudinal spread increases only logarithmically with energy. 

In Appendix A we give parameterizations of the frequency spectrum of the 
Cherenkov electric field emitted in hadronic showers, for practical applications, 
and we discuss their range of validity. The magnitude of the electric field is 
given as a function of proton energy, observation angle and frequency. As an 
example we show in Fig. [10] the electric field spectrum at the Cherenkov angle 
in a 100 PeV proton shower, compared to the corresponding parameterization 
as given in Appendix A. As stated in Appendix A, one can see that the 
parameterizations work very well (accuracy ~ 1%) for frequencies up to the 
frequency z/max? at which the spectrum is maximum for each observation angle. 
For frequencies above i^max, the accuracy worsens gradually and reaches ~ 5% 
at z/ = 2i/max, for observation at the Cherenkov angle 6 = 9c, and ~ 15% at 
angles 6 = 6ci: 10°. In the top panel of Fig [10] we also show parameterizations 
of the frequency spectrum of 100 PeV proton-induced showers as obtained 
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Fig. 8. Average frequency spectrum of the Cherenkov radiation obtained in 
ZHAiRE S simulations of electron (solid blue line) and proton (dashed red line) 
showers with primary energy £^0 = 1 PeV (top panel), Eq = 100 PeV (middle 
panel) and Eq = 10 EeV (bottom panel). The spectrum is shown at four observa- 
tion angles with respect to the shower axis. The RMS of the 20 simulated showers 
is also shown. 
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Fig. 9. Same as Fig. [HI but showing only proton-induced showers at energies Eq =100 
TeV(long-dashed black line), 10 PeV (dashed red line) and 1 EeV (solid blue line). 
The electric field was divided by shower energy. 

in [33]. Those parameterizations are based on l-dimensional simulations and 
the lateral distribution was only acccounted for in a very approximate way, 
which explains the large discrepancies with the results of this work at the 
Cherenkov angle. However, the fits are in fairly good agreement for angles 
away from the Cherenkov angle since the features of the spectrum at those 
angles are determined by the longitudinal development of the shower which is 
correctly accounted for in [33] as well as in this work. 



3.3 Contribution of charged pions, muons and protons to excess track lengths. 



In a hadronic shower, there are more particles besides electrons that can be 
expected to contribute to the excess charge. In fact, it has been experimentally 
determined [57] in atmospheric showers that cascades develop an excess of 
positive muons over negative muons, stemming from an excess of positively 
charged pions and kaons [SS] • As a result, we expect pions and muons to 
contribute with an excess of positive track lengths to the Cherenkov radio 
emission. 

In Fig. [11] we show the average longitudinal profile of electrons, positrons, 
protons, anti-protons, vr"*", vr", /i"*", /i~, K^ and K~ as obtained in ZHAire S 
simulations of 20 proton- induced showers of ii^o = 1 PeV in ice. It becomes 
apparent that the main contributions to the total projected tracklength, which 
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Fig. 10. Top panel: Frequency spectrum of the Cherenkov radiation obtained in 
ZHAiRE S simulations of a proton shower with primary energy Eq = 100 PeV 
(dashed blue line) at the Cherenkov angle and 5° and 10° away from it. The numer- 
ical results are compared to the parameterisations for proton-induced showers given 
in Appendix A (solid red lines) and in [33] (dotted magenta line). Bottom panel: Fre- 
quency spectrum obtained in ZHAire S simulations of a neutrino-induced shower 
in a neutral current interaction (or a charged current interaction of a muon neu- 
trino) in which the secondaries from the fragmentation of the nucleon carry an 
energy Ej^t = 2 EeV (dashed lines). The frequency spectrum is compared to the pa- 
rameteresation of a proton shower with energy Eq = 2 EeV (solid line) at different 
observation angles. See Section S] for further details. 
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Fig. 11. Longitudinal development of of electrons, positrons, protons, anti-protons, 
vr"'', vr^, fi^, /i~, K^ and K~ as obtained in ZHAire S simulations of proton- in- 
duced showers in ice at 1 PeV. 

is approximately proportional to the area under the longitudinal profile, come 
from charged pions, protons and muons. We neglect the contribution from 
charged kaons, which amounts to ~ 1% of that due to pions. Protons con- 
tribute most to the excess positive tracklength because almost no anti-protons 
are produced in the shower to compensate for them, as can be seen in Fig. [TTl 

In order to quantify the influence of charged pions, muons and protons on 
the track lengths, we computed both, the total and projected track lengths of 
those particles in ZHAire S simulations. As expected, we found that there is 
an excess of track lengths due to positively charged pions, muons and protons. 
For example, in a 10 PeV proton shower: 



Ri 



(Te--Te+) 



0.015 



(3) 



where T denotes the projected track length due to different particles. This 
means that the excess of positively charged pions, muons and protons dimin- 
ishes the excess projected track length of the shower by ~ 1.5% because the 
excess is of opposite sign, compared to the excess of electrons. This contribu- 
tion decreases in magnitude with shower energy in the same fashion as the 
EM energy ratio increases (see Fig. [7| and Table [2]). Also, for a fixed shower 
energy, there is a dependence on the hadronic interaction model used in the 
simulations, as will be discussed in the next section. 
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In the top half of table [2] we show the total and projected track length divided 
by the primary energy Eq due to different particles in proton-induced showers 
oi Eq =1 TeV, 1 PeV and 1 EeV. Note that the e^ track lengths normalized 
to Eq increase with energy, since the energy feeding the electromagnetic com- 
ponent of the shower also increases with the energy of the primary proton. 
This is illustrated by the EM energy fractions shown in the bottom of table [2j 
As discussed before, this increase in the EM energy ratio is what causes the 
proton-induced shower projected excess e^ tracks to increase with energy, as 
can be seen in figure [3 Furthermore, both the total e^ track/£'o and the pro- 
jected e^ track/ii^o are proportional to the EM energy fraction of the shower. 
The ratio between the total e^ track length normalized to Eq, and the EM 
energy fraction of the shower is constant, i.e. {J2Ti.± / Eq) / EM^^"''^^'^°^ ~ 57 
m/TeV. The equivalent ratio for the projected e^ tracks varies only slightly 
between 41 and 45 m/TeV for all the energy range shown in table [H 



Total Track/^o [m/TeV] 


Proj. Track/£;o [m/TeV] 




1 TeV 


1 PeV 


lEeV 


1 TeV 


1 PeV 


1 EeV 


e+ 


1271(145) 


1679(43) 


1747(16) 


1024(158) 


1466(47) 


1539(17) 


e~ 


2629(287) 


3440(86) 


3573(32) 


1801(262) 


2544(79) 


2664(30) 


7r+ 


280(69) 


93(20) 


62(7) 


195(50) 


65(14) 


43.1(5.2) 


7r~ 


277(74) 


90(20) 


60(7) 


184(51) 


59(13) 


39.6(4.8) 


/^+ 


14.7(7.3) 


5.1(1.2) 


3.6(0.6) 


11.0(6.8) 


3.8(0.9) 


2.4(0.4) 


/^" 


14.4(7.2) 


5.0(1.5) 


3.3(0.4) 


9.6(6.1) 


3.6(1.3) 


2.16(0.35) 


V 


85(24) 


27.9(6.3) 


18.6(2.2) 


68(20) 


21.6(4.9) 


14.4(1.7) 


P 


0.4(0.9) 


0.3(0.2) 


0.19(0.03) 


0.3(0.9) 


0.3(0.2) 


0.18(0.03) 





^^fraction protoii-induced 


^^fraction electroii-induced 


1 TeV 


68.7% 


99.6% 


1 PeV 


89.2% 


99.2% 


1 EeV 


92.6% 


98.0% 



Table 2 

Top table: Average total track lengths (top left) and projected track lengths (top 
right) divided by £"0 for different particle species obtained from simulations of 
proton- induced showers of Sq =1 TeV, 1 PeV and 1 EeV in ice, using ZHAire S. 
Also shown in parenthesis are the RMS values for the 20 showers. Bottom table: 
Fraction of shower energy feeding the EM component for proton- and electron- in- 
duced showers of Eq =1 TeV, 1 PeV and 1 EeV. The cuts used in these simulations 
were 80 keV for e and 500 keV for hadrons. 

The very small contribution of pions, muons and protons to the field can be 
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Fig. 12. Average frequency spectrum of the Cherenkov radiation obtained from 
ZHAiRE S simulations of 10 PeV proton-induced showers accounting (red dashed 
lines) and not accounting (blue solid line) for pions, muons and protons in the 
calculation of the electric field. The spectrum is shown at four observation angles 
with respect to shower axis. The RMS of the 20 simulated showers is also shown. 

seen in Fig. [121 where we show the Fourier-spectrum of the electric field in 
proton showers of 10 PeV with and without accounting for the charged pions, 
muons and protons tracks. Although the actual field emitted by the hadrons 
in the shower is very small, compared to the field emitted by electrons and 
positrons, the hadronic component carries a significant fraction of the shower 
energy, especially at lower energies. The importance of this energy balance 
of the shower can be illustrated by the effects of the EM energy fraction on 
the radio emission, as discussed before. Given our results, following the actual 
tracks of hadrons in the shower turns out not to be essential for the field 
calculation, but accounting for the hadronic energy of the shower is. 



3.4 Hadronic interaction model dependence 



Modeling of hadronic showers above PeV energies relies on the extrapolation 
of hadronic interactions to energies and regions in the parameter space of 
the collisions which have not been probed in terrestrial accelerators [59]. We 
have studied the influence of the hadronic model on the coherent Cherenkov 
radiation frequency spectrum. For that purpose we have simulated a set of 
proton-induced showers with the QGSJETOl model [60], and compared the 
spectrum with that obtained using the SIBYLL 2.1 model [HI], which we have 
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been using by default in all other simulations so faiL^. This comparison is 
shown in Fig. [131 The differences between the normalization of the spectra 
in the coherent region obtained with SIBYLL 2.1 and QGSJETOl are small, 
with the QGSJETOl spectrum smaller than that predicted with SIBYLL 2.1 
by up to ~ 4% around the cutoff frequency at the Cherenkov angle. 

It is well known that QGSJETOl produces more muons and pions per inter- 
action than SIBYLL 2.1 ^63j and also protons. For instance, we have obtained 
the longitudinal profile of the number of pions in proton-induced showers, 
and observed that the number of pions at the maximum is ~ 24% larger in 
QGSJETOl than in SIBYLL 2.1 p3]- Also more protons, pions, and muons 
per interaction means less electrons and positrons in the shower, decreasing 
the negative excess track (denominator of Eq. ([3]) ). As a consequence, we ex- 
pect a difference in the ratio Rt predicted by both models. While in SIBYLL 
2.1 i?r ~ 0.012 at 1 EeV, QGSJETOl predicts a larger ratio Rt ~ 0.014, 
the increase being dominated by the proton projected track-length. Also as a 
consequence, the EM energy ratio is larger (93%) in SIBYLL when compared 
to QGSJET (91%) at 1 EeV. This is consistent with the fact that the electric 
field in the fully coherent region is slightly smaller when using QGSJETOl 
than when SIBYLL 2.1 is adopted. We remark here that the smallness of Rt 
reflects that the tracks of protons, pions and muons are mostly irrelevant to 
the field calculation, and this conclusion holds regardless of the details of the 
high-energy hadronic interaction model. In fact, we have also studied the effect 
of the low energy hadronic models on the radio emission, and our conclusions 
remain the same. 

There are also differences in the frequency at which the electric field intensity 
is maximum at each observation angle away from the Cherenkov angle. These 
stem from the fact that the longitudinal profile of the excess negative charge 
predicted by SIBYLL 2.1 is longer than in QGSJETOl, leading to higher 
frequencies at the maximum when using QGSJETOl. Also, the lateral profile 
is slightly flatter in QGSJET, leading to a slightly smaller cut-off frequency 
at the Cherenkov angle. However, these differences are small. 



4 Summary and conclusions 



We have presented ZHAire S, a Monte Carlo code that merges the high energy 
hadronic interaction and tracking capabilities of AIRES ^44], and the dense me- 
dia propagation capabilities of TIERRAS [45j with the precise low energy e^ 
tracking and radiation calculation capabilities of ZHS pi25] . This combination 



^ We use the SIBYLL 2.1 and QGSJETOl models since recent data favors them 
over other models 1621. 
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Fig. 13. Average frequency spectrum of the Cherenkov radiation obtained from 
ZHAiRE S simulations of 1 EeV proton-induced showers using QGSJETOl (red 
dashed lines) and SIBYLL 2.1 (blue solid line). The spectrum is shown at four 
observation angles with respect to shower axis. The RMS of the 20 simulated show- 
ers is also shown. 

allows a precise full simulation of radio emission of electromagnetic, and for 
the first time also hadronic showers up to EeV energies. ZHAire S has been 
used in conjuction with the TIERRAS package [IS] to simulate showers in 
dense media, obtaining the radiation emitted due to the Askar'yan effect. 

We have compared ZHAire S results for electromagnetic showers in ice against 
ZHS up to 10 EeV, obtaining a good agreement between them. In the case of 
hadronic showers, we have compared the results from ZHAire S and GEANT4 
up to 90 TeV [28], the highest energy for full simulations of hadronic show- 
ers with GEANT in the literature, also with very good quantitative agree- 
ment. Unlike GEANT4, ZHAire S is capable of simulating the emission from 
hadronic showers up to EeV energies. 

By comparing the results of thinned ZHAire S simulations with un-thinned 
ones up to 1 PeV, we obtained thinning parameters which give a good com- 
promise between accuracy and CPU time. We then tested the applicability of 
these parameters at higher energies, by comparing our results against ZHS, 
which has its own, well established [32], thinning algorithms and parameters. 

We confirmed the expectation that the EM energy ratio of hadronic showers 
keeps approaching that of EM showers at the highest energies, increasing the 
magnitude of the electric field (normalized by £"0) as energy grows. We also 
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showed that, in the case of UHE purely EM showers, the EM energy ratio de- 
viates from a hnear dependence on energy, due to photohadronic interactions. 

By comparing EM and hadronic showers up to 10 EeV we have found that, 
at the Cherenkov angle, the cut-off frequency tends to be smaller in hadronic 
showers and increases slowly as the energy rises because the lateral distribution 
becomes narrower. For angles away from the Cherenkov angle, as expected, 
the cut-off frequency decreases with energy, due to the logarithmic growth 
of the longitudinal profile. In contrast, away from the Cherenkov angle, the 
cut-off frequency for EM showers decreases rapidly above PeV energies due 
to the LPM effect, which is much more pronounced in EM showers than in 
hadronic ones. We also found that the fluctuations in hadronic showers are 
almost independent of energy, while in EM showers they grow rapidly, as the 
LPM effect becomes important. 

We have also analyzed the influence of charged pions, muons and protons on 
the radio emission of showers. We found that most of their contribution is due 
to the excess of protons, with pions being the next contribution in importance, 
that induce a decrease of the (total) negative excess track length of the shower 
of ~ 1 — 2% above 1 PeV. Although the actual field emitted by the hadrons in 
the shower is very small, the hadronic component carries a significant fraction 
of the shower energy. The importance of this energy balance is illustrated 
by the effects of the EM energy fraction on the radio emission. Given our 
results, following the actual tracks of hadrons in the shower turns out not to 
be essential for the field calculation, but accounting for the hadronic energy 
of the shower is. 

We have compared the results obtained using SIBYLL 2.1 [61j and QGSJETOl 
|60] . and found that there is a small dependence of the normalization of the 
emitted electric field on the hadronic model, up to ~ 4% lower in QGSJETOl 
at the Cherenkov angle. This is due to the higher number of protons, pions 
and muons when QGSJETOl is used, which increases the positive projected 
excess track-length of protons, pions and muons, with protons dominating the 
difference. Also more protons, pions and muons means means less electrons 
and positrons in the shower, decreasing the negative excess track. The net 
effect is a lower negative projected excess track as a whole, diminishing the 
field. At angles away from the Cherenkov angle, the longitudinal excess profile 
becomes important, and the shorter showers obtained using QGSJET lead 
to the maximum of the emission at higher frequencies, when compared to 
SIBYLL showers. 

Finally, we have used our simulations of the emission due to the Askar'yan ef- 
fect in hadronic showers to obtain new parameterizations of the radio pulse fre- 
quency spectrum which are described in Appendix A. We note that our results 
for proton showers can also be applied to approximately model the Cherenkov 
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emission from hadronic showers induced in high-energy neutrino interactions. 
This is shown in the bottom panel of Fig. [10] where we have plotted the 
frequency spectrum obtained in ZHAire S simulations of a neutrino-induced 
shower in a neutral current interaction (or a charged current interaction of a 
muon neutrino) in which the secondaries from the fragmentation of the nu- 
cleon carry 2 EeV of energy in total. The products of the neutrino-nucleon 
interaction, mainly pions, were obtained with the HERWIG interaction Monte 
Carlo code [M]- The frequency spectrum obtained in ZHAire S simulations 
is compared to the parameterization of a proton shower with energy Eq = 2 
EeV as given in Appendix A. 
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6 Appendix A 



In this Appendix we give parameterizations of the magnitude of the Cherenkov 
electric field in proton-induced showers in ice for practical appications, as 
obtained in simulations with ZHAire S. 

This is the first time parameterizations of the electric field in hadronic showers 
are presented. The parameterizations given here complement those presented 
in [32] for purely electromagnetic showers. However in [32] photoproduction 
interactions were not taken into account, and for this reason the parameter- 
izations presented there for electromagnetic showers should be less accurate 
for E > 10^^ eV. This important subject will be further studied in a future 
work. 

The physical basis for our parameterization of the radiated spectrum is the 
"box model" of shower development [25|lT|rr|l^ . In this model, the distri- 
bution of particle tracks has a characteristic length L, proportional to the 
radiation length Xq, and a width R, proportional to the Moliere radius Rm- 
The model motivated a functional form for the radiated spectrum -E(z/) in 
electromagnetic showers, given by [32] : 



r\EiEo,e,u)\ = AiEo,e,u) x dLiEo,9,iy) x dRiEo,e,u) (4) 

where A is a coherent amplitude that increases linearly with frequency, which 
is multiplied by lateral and longitudinal decoherence factors dn and dL, defined 
in such a way so that d^^L < 1 and dn^L — >■ 1 as i/ — )■ 0. Also, as 6* — )■ 9c, 
where 6c is the Cherenkov angle, c^l — )■ 1, since interference over the width of 
the shower dominates. In general. A, d^ and dji depend on the shower energy 
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Eq, frequency u, and observation angle 6 with respect to the shower axis, r is 
the (far-field) observation distance. 

In this paper we adopt this same functional form for the radiated spectra in 
proton-induced showers. For the decoherence factors we use the same func- 
tional forms as in [7], namely: 



where un and ul are characteristic frequencies at which lateral and longitu- 
dinal decoherence become important, and a and f3 give the strength of the 
decoherence. The frequencies ul and z/^j are inversely proportional to L and 
R, respectively. The amplitude A of the fully coherent (low-frequency) com- 
ponent is proportional to the total excess track length T, which in hadronic 
showers is known to deviate from a linear behavior with shower energy (see 
Fig. [7]). The length L, width R, and excess track length T of the shower can 
be related to properties of the interaction medium [25|41II31|I7] : 



urie) ^ ^ = ^ (6) 

^' L\l~ncose\ kLXo\l-ncos9\' ^' 

MO = 0c) - ^ = /. ," , (T) 

R krRm V^ - 1 

A{Eo,e,iy) = kEiyTsinO ^ fcg „° ° t/MHzSin6' , (8) 

Ec P 

where j^mhz is the frequency in MHz, fc^ has units of V cm~^ MHz~^ and c = 
3 • 10^° cm s~^ is the speed of light. Xq = 36.08 g cm~^, Rm = 10.57 g cm~^, 
n = 1.78, p = 0.924 g cm~^ and Ec = 73.1 MeV are respectively the ra- 
diation length, Moliere radius, refractive index, density, and critical energy 
of ice. Throughout, we use the approximation ur^O) = i'r^Oc), since lateral 
decoherence is only important near the Cherenkov angle. 

The quantity Jce is also energy dependent, since the track length due to the 
excess negative charge deviates from linearity in hadronic showers (see Fig. [7]). 
To fit kE^Eo) we use: 



where kEfl, kE,i and Ee are constants given in Table [31 

Since L is dependent upon the shower energy Eq, ki is in fact also energy- 
dependent, although with a much weaker dependence than in electromagnetic 
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kEfl [V cm-1 MHz-2] kE,i logio(^i=;/eV) 

4.13 lO"^*' 2.54 10.60 

Table 3 

Fitted proton-induced shower parameters in ice, as defined by Eq. ([9]). 

showers, because hadronic showers are not so strongly affected by the LPM 
effect. For this reason we expect k^ to be well fit in the whole energy range 
by the relation: 

kL = kL,o [-^) , (10) 



where ki^, E^ and 7 are constants given in Table IH 

We also found a weak dependence of kpi with energy, which we further param- 
eterize as: 

kR = kRfi + tWi\i\ I , (11) 



where kn^Q, kpt^i and Er are constants given in Table El 

Finally, a and /3 are found to be practically independent of shower energy, 
and they are given in Table [61 

Our fitting procedure is the same as in [32]. We reproduce it here for self- 
consistency of the paper. For each shower of energy Eq, we start by obtaining 
A{6, v) and (Ir from a fit to the spectrum at the Cherenkov angle 6c where 
dL = 1. We then obtain Ize by substituting A{6, v) in Eq. ([H]). Also from the 
fitted (Ir in Eq.(l5]), we obtain parameters pr and a, and relate vr to k-R via 
Eq. (I7j). Fixing A and (Ir to their fitted values, we then allow d^ to vary on fits 
of Eq. (jl]) to the simulated field at various angles away from the Cherenkov 
angle, obtaining dL{9, v) for each angle. From these we can obtain viid) and /3 
from Eq. ([5]), and ki from Eq. ([6]), respectively. We repeat this process many 
times, varying the energy Eq in multiples of 10 over the energy range from 1 
TeV to 10 EeV. By using the values of kE^EQ), ki^EQ) and Jcr^Eq) from fits at 
these different energies, we obtain the parameters of Eqs. ([9]), fITOl) and flTTj) . 
that define the energy dependence of kE, ki and kR, and are given in Tables 
[31 [H and [5l respectively. The fitted parameters a and (3, practically energy 
independent, are given in Table [B 

While all the fitted parameters vary on a shower-to-shower basis, variations 
tended to be small, so that approximating kE ~ ^e, f3 ^ (3, etc., is appropriate 
for individual showers. This is in contrast to previous parameterizations of 
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Table 4 



khfi 


7 


logio (i^L/eV) 


31.25 


3.01 10-2 


15.00 


3ed shower 


parameters in ice, as defined by I 


fciJ.o 


kR,i 


logio(i^ij/eV) 


2.73 


1.72 


12.92 



Table 5 

Fitted proton-induced shower parameters in ice, as defined by Eq. ([TT 



a P 



1.27 2.57 

Table 6 

Fitted proton-induced shower parameters in ice, as defined by Eq. (0). 

electromagnetic showers, in which k^ varied strongly from the mean fitted 
values Jcl, due to variations in the longitudinal spread of the showers caused 
by the LPM effect. This is not the case in hadronic showers and we take 
ki ~ kL. 

Finally, the fits have an accuracy of ~ 1% for frequencies up to the frequency 
i^max; ^t whlch the spectrum is maximum for each observation angle. For fre- 
quencies above t'max, the accuracy worsens gradually and reaches ~ 5% at 
u = 2z/inax; ^OT observatlou at the Cherenkov angle 6 = Oq, and ~ 15% at 
angles 6 = 60 ±10°. 
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